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Abstract — In this work, we consider the effects of gas slippage and wellbore storage in shale gas 
reservoirs. We use the Finite Difference Method for the discretization of the nonlinear governing equa¬ 
tion, and the iterative Gauss-Seidel method is applied to obtain the solution of the algebraic system. 
We also perform pressure tests in the producing well through the use of numerical simulation using 
cylindrical coordinates. The results, obtained in the well testing analysis context, show the relevance of 
the introduction of the slip, formation damage and wellbore storage on the flow simulation in shale gas 
reservoirs. 
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I. INTRODUCTION 

Archaeological records show that ancient civiliza¬ 
tions in different parts of the world used oil. Incas and 
Babylonians used the oil to lay bricks, to pave roads, 
and to waterproof ceramic artifacts [6, 34], However, 
it was only in the 20th century that oil began to stand 
out as economic activity on the world stage, with the 
emergence of the large automotive industries and oil 
companies, with the latter beginning the search for 
more oil reserves. On the other hand, natural gas 
has been in the background for a long time, mainly 
due to its specificities concerning transport and stor¬ 
age. However, this situation has changed, mostly in 
the transition between the 20th and 21st centuries. 

The market for natural gas of fossil origin has 
grown worldwide for several reasons, such as the dis¬ 
covery of new deposits, the lower degree of pollution 
that its burning causes, when compared to oil, and 
its diversity of applications. The economic relevance 
of natural gas was evident in the years 1973 and 
1979 when the Organization of Petroleum Exporting 
Countries (OPEC) raised oil prices to a level that led 
to changes in hydrocarbon exploration and produc¬ 
tion and the general dynamics of generation and con¬ 


sumption of energy in the world. As a consequence, 
activities associated with natural gas increased [4], In 
this context, Fig. 1 presents the Brazilian energy ma¬ 
trix in the 1970s and 2010. In this figure, we can see 
that the consumption of natural gas increases from 3 
to 9%. 

In Brazil, in 1941, natural gas began to appear in 
the economic scenario due to the discovery of the 
first commercial oil field in Candeias, Bahia, the in¬ 
dustrial sector being its first destination. Due to the 
adverse effects of the two oil shocks, there was an 
increase in the production of natural gas due to the 
exploration and development of the Campos Basin 
in Rio de Janeiro [9]. At the end of the 1990s, with 
the completion of the construction of the Bolivia-Brazil 
gas pipeline (GASBOL) [7], the import of natural gas 
from Bolivia began, and it consolidated itself as a rep¬ 
resentative energy source in the Brazilian energy ma¬ 
trix, leaving aside the oil by-product label. 

Currently, with the discoveries of pre-salt reserves 
in the Santos Basin [32], interest in natural gas re¬ 
serves in Brazil has been growing even more rapidly 
in recent years. Figure 2 shows the growth of gas pro¬ 
duction over the years and the dependence on Bo- 
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livian gas imports, which is necessary to supply na¬ 
tional consumption. However, the Energy Research 
Company (EPE) disclosed, through its National En¬ 
ergy Plan [14], that the national natural gas supply will 
be adequate to meet the demand of all the commit¬ 
ments assumed until 2030, thus making self-sufficient 
Brazil in the sector. 
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Fig. 1: Brazilian energy matrix in the 1970s and 
2010 [15]. 
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Fig. 2: Gas production and consumption in Brazil 
between 2008 and 2017 based on the 
National Energy Balance [ 19], 


1.1. Obtaining natural gas 

We define a reservoir of natural gas as a de¬ 
posit that presents a mixture of hydrocarbons in the 
gaseous state [17], As stated by Ezekwe [17], we can 
classify gas reservoirs as a reservoir of dry gas, wet 
gas, or retrograde condensation gas, according to the 
type and behavior of the fluid. When there are no 
molecules heavy enough to form liquid hydrocarbons 
after surface separation processes, we say that the 
reservoir is to be dry gas. We call the deposit a wet 
gas reservoir when it produces a certain amount of liq¬ 
uid in the well and/or in the surface facilities. In turn, 
as we produce gas from a retrograde condensation 
gas reservoir, the pressure decreases, and the tem¬ 
perature remains approximately constant. Thus, there 
is a liquid phase formation in the reservoir through the 
condensation of the gas. 

Another classification that can be applied is that of 
conventional and unconventional reservoirs. In con¬ 
ventional gas reservoirs, it is possible to produce from 
vertical wells while it is required specific production 
techniques for unconventional reservoirs. For exam¬ 
ple, we can cite the horizontal wells and hydraulic frac¬ 
turing since there is higher resistance to the gas flow 
due to the lower absolute permeability of the reser¬ 
voir. However, with the advancement of technological 
development, the gas of unconventional origin began 
to gain prominence in the world economic scenario. 
India [21], Brazil [27], and USA [40] are among the ten 
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nations with the largest reservoirs of shale gas in the 
world, a type of unconventional reservoir. The trend 
is that gas of unconventional origin will be essential in 
the world economy in the coming decades. 

Unconventional gas reservoirs can be of different 
types: deep gas reservoirs (located beyond 4,500 
meters in depth); low permeability (tight gas); shale 
gas; gas adhered to coal veins (coalbed methane); 
gas from pressurized zones (at very high pressure 
when compared to other reservoirs with the same 
depth); and underwater and Arctic hydrates (methane 
hydrates) [2]. The focus of this work, however, is the 
shale gas reservoirs. 

1 . 2 . Shale gas reservoirs 

The shale gas reservoir sedimentary rock has 
a fine granulometry and very low permeability [22], 
Also, the shales have quite variable mineralogical for¬ 
mation, with a predominance of brittle minerals such 
as quartz, carbonates, and feldspars [20]. Although 
it has characteristics of unconventional reservoirs, its 
exploitation has increased over the years. 

Hasan et al. [21] has shown that two-thirds of the 
world’s hydrocarbon reserves are unconventional and 
that this fact is directly related to the growing impor¬ 
tance of shale gas in the world energy matrix. On 
the other hand, Gomes [19] claims that technological 
advances in the sector are allowing its production to 
become attractively economical. It has been drilled 
more than 50,000 wells over two years in gas reser¬ 
voirs in the United States [28], making the country’s 
hydrocarbon production almost double. 

Countries such as the USA and Brazil have a large 
number of gas reservoirs, which is very important in 
economic terms and of influence in the international 
socio-political framework. As a result of the growth 
in non-conventional gas exploration, these countries’ 
dependence on the world’s largest oil and gas produc¬ 
ers, such as Venezuela and Russia, can end. 

Despite the environmental concern with fossil fu¬ 
els and the rapid growth in the use of fuels from re¬ 
newable sources, Jia et al. [22] believe that fossil en¬ 
ergy should still account for 78% of global energy 
consumption in the year 2040. Even though natural 
gas of fossil origin is not renewable, technological ad¬ 
vances in the area and the growth of discoveries of 
shale gas reservoirs guarantee its use for still many 
years. Knudsen et al. [24] show that in recent years 


the United States has been reducing the use of coal- 
fired power plants and replacing them with natural gas 
in the generation of electricity. In part, we can explain 
this by the desire to reduce emissions of polluting car¬ 
bon dioxide, which harms the environment so much. 

Wellbore Testing Analysis studies the pressure 
and flow changes as a function of production time, 
through measurements at the bottom of the well and 
the flow at the surface. From the measured pres¬ 
sure response, it is possible to determine the reser¬ 
voir properties useful for production planning [41], In 
a well test, a transient pressure response occurs due 
to a production/injection flow. Depending on the ob¬ 
jectives of the test, we record the response of the well 
over a relatively short period when compared to the 
productive life of the reservoir. In this work, in addi¬ 
tion to a sensitivity analysis, well pressure tests using 
numerical simulation is also done. 

II. POROUS MEDIA GAS FLOW 

Petroleum is a mixture of hydrocarbons, which has 
its physical state-properties determined, in general, 
by its composition, temperature (T), and pressure (p). 
According to Ezekwe [17], oil is the part that remains 
in the liquid state when a mixture of hydrocarbons 
is brought from the reservoir conditions to the sur¬ 
face conditions, while the natural gas presents the 
gaseous state in the surface conditions. Under reser¬ 
voir conditions, natural gas can be present in gaseous 
form or dissolved in oil. 

Regarding the composition of natural gas, the 
amount of each component can vary depending on 
the type of reservoir and its characteristics, for exam¬ 
ple, the location (land or sea), the type of soil, and 
the geological formation process of the basin, among 
other factors [30]. However, as can be seen in Ta¬ 
ble 1 which presents the typical composition of a nat¬ 
ural gas reservoir, it is evident that the primary com¬ 
ponent is methane, which may represent an amount 
of 70 to 98% of the total natural gas, and in smaller 
quantities, considered as impurities, carbon dioxide, 
hydrogen sulfide, and nitrogen. Besides temperature 
and pressure conditions, another fundamental param¬ 
eter for calculating the properties of the gas is its rel¬ 
ative density (or specific gravity), 7, which is the ratio 
between the molecular mass of the gas, M, and the 
molecular mass of air, M air . In this work, it is con¬ 
sidered a reservoir of dry gas, produced without the 
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appearance of liquid at any time of production. 

Table 1: Typical chemical composition of natural gas 


Component 

Composition 

n 2 

up to 15% 

C0 2 

up to 5% 

h 2 s 

up to 3% 

He 

up to 5% 

ch 4 

70-80% 

c 2 h 6 

1-10% 

C3H8 

up to 5% 

C4H10 

up to 2% 

CsH 12 

up to 1 % 

CeHi 4 

up to 0.5% 

C7+ 

up to 0.5% 


For gas mixtures of hydrocarbons, we use the 
pseudo-critical pressure and temperature coordi¬ 
nates, p pc and T pc , respectively, to determine the so- 
called pseudo-reduced coordinates [26], and we use 
them to calculate the physical properties of natural 
gas. 

Sutton [38] presents, depending on the gas den¬ 
sity, the correlations that we apply here to obtain the 
pseudo-critical pressure and temperature, which are 
fundamental for the reservoir simulations [16]. We use 
them in determining, for example, the compressibility 
factor (Z), volume formation factor (B) and viscosity 
(p) [37], From Z and the universal gas constant, R, it 
is possible to determine the density, p, from the equa¬ 
tion of state for a real gas, p = pM/ZRT. 

The gas volume formation factor is the relationship 
between the volumes it occupies under reservoir con¬ 
ditions (V) and standard conditions (V sc ) (pressure, 
p sc , and temperature, T sc , in standard conditions) [16]. 
Thus, B = p sc ZT/pT sc whereas Z ac « 1. 

On the other hand, we calculate the viscosity of 
natural gas using the correlation suggested by Lee 
et al. [25], widely used in reservoir simulation. 

Here, the effective porosity (4>) varies depending 
on the pressure [16]: 

(j) = (j) 0 [1 + (p - p°))\ ( 1 ) 

where 4>° and p° are, respectively, the porosity and 
pressure in the reference conditions, is the coef¬ 
ficient of compressibility of the rock, and we assume 


that the compressibility of the rock is small and con¬ 
stant. 

In addition to porosity, the economic viability of 
a reservoir also depends on the permeability of the 
rock. This property is a measure of a porous mate¬ 
rial’s ability to allow fluids to pass through its pores. 
We usually represent the absolute permeability by the 
tensor k. 

2.1. Slip in gas flow 

Studies and predictions about the flow of gas in 
porous media are more difficult to carry out than those 
of liquid because the gas properties generally depend 
more strongly on pressure and also due to the differ¬ 
ent mass transport mechanisms that can be present 
[28, 31 ]. Therefore, in some cases, the classic Darcy’s 
law does not adequately describe the flow physics, 
and experimental data suggest corrections for the cal¬ 
culation of permeability and, thus, we introduce a 
modified Darcy’s law [28] 

v = ^h(Vp-pgVD), ( 2 ) 

h 

where k a is the apparent permeability tensor, v is the 
surface velocity of the fluid, g is the acceleration of 
gravity and D is the depth. 

Specifically for gases, the slip flow regime occurs 
when the average free path of the gas molecules has 
a scale comparable to the pore size [18]. So, both 
the reservoir and fluid properties influence the deter¬ 
mination of apparent permeability. We can mention, 
among the non-Darcy effects, that we can incorporate 
in the apparent permeability, high flow rates (inertial 
and turbulent effects), non-Newtonian fluid flow (for 
liquids), and slip flow, which occurs only for gas under 
certain reservoir conditions [5]. 

When the fluid is a gas, the Klinkenberg effect 
shows that the permeability measurements made in 
the laboratory result in values higher than the abso¬ 
lute values, due to the slipping of the gas on the walls 
of the porous medium. This slip results in a higher 
flow and leads to a correction of the apparent perme¬ 
ability [23], 



where b is the Klinkenberg parameter and k the abso¬ 
lute permeability tensor. 
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In reality, the mass is transported in the porous 
medium by a variety of mechanisms, one of which 
is the so-called Knudsen diffusion. The Knudsen 
number measures the relationship between the mean 
free path of the molecules, A, and the character¬ 
istic pore length, R h , so that Kn = A /R h , where 
A = (, ii/p)y/nZRT/{2M ), R h = 2yfcsjkj^, and r 
is the tortuosity of the porous medium. 

The slip flow regime occurs for 10~ 3 < Kn < 0.1 
and we can also introduce a different model for deter¬ 
mining apparent permeability [28] 

k “ = ( 1 + rr^) k = /( ^ )k ' (4) 


2.2. Governing equation 

In obtaining the partial differential equation (PDE) 
that governs the isothermal flow of a gas in a porous 
medium, we employ the mass conservation equation 
and the modified Darcy’s law. We also take into ac¬ 
count the effects of slippage, wellbore storage and 
formation damage, disregarding the phenomenon of 
gas adsorption, the gravitational force and non-Darcy 
behaviors related to inertial effects. 

As it takes time to the hydrocarbons in the reser¬ 
voir to reach the surface, in the first moments, we 
produce the fluid initially stored in the well. This ef¬ 
fect is called wellbore storage. The formation damage 
concerns the reduction of permeability in the region 
close to the well, caused by wellbore fluids used dur¬ 
ing drilling and completion. 

According to Li et al. [29], mass conservation for 
the flow of gas in porous media can be described, ex¬ 
cluding the adsorption effects and source terms [3], 
by 



account the fluid and rock properties [13] 


d_ (<j>_ 
dt \B 


1 dcj) dp d / 1 A dp 
B dp dt + dp \B) dt 
dp 
dt 


’C0^° 

, d 
+ 4 >— 

'i'll 

|_ B 

v dp 

\B) \ 


= r, 


dp 
'dt ’ 


(7) 


where we also employed Eq. (1). 

To study the flow dynamics in the region close 
to the producing well, we assume a two-dimensional 
flow in cylindrical geometry in the r^-plan and a diag¬ 
onal permeability tensor, so that: 


Id / k or dp A d / k az dp\ dp 

r dr \pBdr) dz \pBdz)~ p dt' 1 ’ 

where k ar and k az are the apparent permeabilities in 
the r- and z- directions, respectively. 

As we are considering the effects of storage in the 
well, the production flow in it is given by [29] 

,, _ s, dp w f 

QSC Qsc “I” C sc “77 1 (9) 

dt 

where q sc is the flow from the porous medium, C sc is 
the storage coefficient (that already incorporates B) 
and p w f the pressure in the well. The flow rate q sc is 
calculated by [35] 


Qsc d w (p Pwf') 5 ( 13 ) 

where J w is the productivity index. 

Finally, as an initial condition we impose 


p{r, z,t = 0) = p ini (r, z) = p ini , 


( 11 ) 


where p ini represents the initial pressure before the 
reservoir undergoes any changes due to fluid produc¬ 
tion/injection. 

On the other hand, the external boundary condi¬ 
tions are of null flow at the external borders, 




( 12 ) 


Then, replacing Eq. (2) in Eq. (5) and ignoring the 
effects of gravity due to the low specific gravity of the 
gas and the thickness of the reservoir, 


where L z is the thickness of the reservoir and r e is 
the outer radius of the reservoir. In turn, for the inter¬ 
nal boundary condition, 


V • 



d_ 

dt. 



( 6 ) 


We can rewrite the term d(4>/B)/dt if we take into 



QscBp 

Q.TTkayh'Pyy 


(13) 


where h is the thickness of the production region con¬ 
sidered and r w is the radius of the well. 
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III. NUMERICAL RESOLUTION 
METHODOLOGY 

We employ a computational mesh of centered 
blocks [1,5,10,16] and the cylindrical coordinate sys¬ 
tem (r and z), Fig. 3, and we also assume the angu¬ 
lar symmetry of the flow. We obtain the numerical 
solution in the nodes of the computational mesh, lo¬ 
cated in the centers of the cells, with n r and n z being 
the numbers of cells in the r- and z- directions, re¬ 
spectively. We also use fractional indexes i ± 1/2 and 
k ±1/2 to indicate cell interfaces of the computational 
mesh. 



i.k -1 

i+l.A-l 

i+2.A-l 

- 

- 

r+l.* 

1+2.* 


i+l.JN-1 

/+2.A+1 



Fig. 3: Two-dimensional cylindrical domain 
discretized. 


So, for the cell i, k and the time n + lwe can write 


that 


, dp\ 
p dt) 


n+1 


i,k 


= v ilk 

d_ 

+ dz 


1 d 


dp 


[ T r + | Ar 

r dr 
dp 


T 


dz 


dr , 

- n+1 

Az 

- i,k 


(14) 


where the governing equation was multiplied by the 
volume of the cell, 1+, and T stands for the transmis- 
sibility. 

Following the techniques traditionally applied in 
reservoir simulation [5, 16] and employing centered 
differences 

d 


Vi k q 
dr 


-c 


Ar. 


i k 

T ^ 

R r <~v 

or 


T, 


dp 

dr 


and 


V 


i,k 


S- frS) 


dz 


n+1 


Vi 


i,k 


dz J i 


T 


A z i<k 

dp 


' dz 


i~b,k 


1+ 

oz 


i.k— i 


i+^,k 

n+1 


(15) 


i,k-\ 

n+1 


(16) 


Similarly, for the discretization of pressure deriva¬ 
tives in the r-direction, 


/ dp N 

n+1 

) 

T) n+1 - T) n+l 

Pi+l,k Pi,k 

Ur, 

*+|,fc 

^ r i+^,k 

and 



/ <9+ 

n+1 

) 


Ur, 

i-\,k 

i—i|CS 

•i 

<1 

1 


(17) 


(18) 


where A r i± y 2 ,k is the distance between nodes of 
cells i, k and i± l,k. We can obtain approximations 
for the derivatives in the z-direction analogously. 

In the case of the accumulation term, we employ a 
conservative expansion [16] 


T pt,l = Vi,k 


1 d(j> 
B n dp 




(19) 


where At is the time step. 

Next, we obtain a totally implicit formulation for 
Eq. (14) using a backward Euler approximation, 


(dp\ n+1 

\ dt )i,k At 


( 20 ) 


3.1. Grid refinement 

For an accurate determination of the pressure gra¬ 
dient, we use to construct the non-uniform mesh in the 
r-direction (where we suppress k for simplicity of no¬ 
tation) [16], 



such as: 


1. we space pressure calculation points using 


T’i+l — Ot-lgTi 

where i = 1, 2, ...,n r - 1; 

2. we define cell boundaries through 

_ r i+ 1 - n 
r * +1/2 log e (n+i/n) 

where i = 1, 2 ,..., n r - 1; and 

3. we calculate cell volumes employing 

r 2 _ r 2 
2 _ 1 i +1 '» 

" I+1/2 “ In (rf + 1 /r?) ‘ 


( 22 ) 


(23) 


(24) 
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In the case of the ^-direction, we use Az = L z /n z . 
We now introduce the definitions of transmissibili- 
ties: 


r+s, 

r,i± ^ 


r,i±-% ,k 

fiB 


n+1 

z±^,/e 


(25) 


and we chose the method of Picard [33]. So this re¬ 
sults in 


T, 


v,n -\-1 
i,k- 1/2 


-u+l.n+1 

Pi,k -1 


+ TV 


v,n+l 
i—1/2,/c 


-u+l.n+l 
Pi-i,k 


T r 


T?,n+1 

i+l/2,fc 


u+l,n+l 

ft+l,* 


+ TV 


i/,n+l 

i,/c+l/2 


u+l,n+l 

Pi.fc+1 


and 


T ”tLi = 


'G n 1+, 

2:,z,/c±^ 

pB 


n+1 


i,/c± i 


(26) 


where, considering that the total angle is 2n and the 
apparent permeabilities [16], 


G n+ + , = 

,k 


X CLTi,k 


log e 


nAz k 


n i+h 


Wi + I,fc 


log e 


n+\ \ 

ri +i) 

(27) 


and 


G"+! , a 

Z,l,K ±^ 


( r h- r U) 


~fc+i - Zk Z k+ 1 -Z k+ x' 


(28) 


x a.Zi, k 


az i,k -\-1 


and we apply arithmetic mean to calculate the fluid 
properties [16]. 


- T 

+ 


v,n -\-1 
i,k- 1/2 


+ T r 


u,n+1 

+ 


i—1/2 ,k 

VAf/ 


v,n -\-1 

v,n -\-1 

+ 


i+l/2,k 

i,fc+l/2 

/ p \ v,n+l 
' 1 P \ „ 

n 


i,k 


v,n+l 

i,k 


(32) 


where the index w refers to the previous iterative level 
when obtaining the pressure. As we can see, the 
properties and coefficients of the discretized equation 
are determined at the iterative level v and used after 
to calculate the new pressure values at the iterative 
level v + l,n + 1. 

In the case of cells in direct contact with the well, 
we determine its pressure (p w f) through Eq. (10) for a 
prescribed flow rate. 

In solving the linearized system, to obtain the pres¬ 
sures in the porous medium and the well, it was cho¬ 
sen to use the Gauss-Seidel iterative method [16, 36]. 


3.2. Numerical approximation for storage 

When we consider the wellbore storage effect, the 
total flow rate is determined as proposed by Li et al. 
[29] and Tavares [39]: 


Q sc — 

with 

Jw - - 


i=K 2 /„"+!_ tfi ' 


k A T 

r ' j ar L - x z- L sc 


PscT In 


(29) 


(30) 


where 


r 0 = ^/r 1+ ir w . 


(31) 


3.3. The solution of the system of equations 

After the process of discretizing the governing 
equation, as the resulting algebraic equations are 
non-linear, we must apply a linearization technique, 


IV. NUMERICAL RESULTS 

In all numerical simulations, we adopted an ar¬ 
rangement consisting of a vertical producer well, of 
length L wf and centered on the r0-plane, a maximum 
production time {t max ), and an initial time step (A t ini ). 
The time step can vary depending on the growth rate 
(5At)- This methodology allows a growing time step, 
and we use it until the final time step (A t max ), pre- 
established, is reached. It is of general use in reser¬ 
voir simulation [1, 16]. We aim to enhance the accu¬ 
racy of the calculated well pressure in the initial mo¬ 
ments when the pressure drop is more accentuated. 

We obtain the results using a standard set of data, 
based on the non-Darcy model discussed in Li et al. 
[29], in which the authors incorporated the slippage 
effect into the apparent permeability. Table 2 shows 
the parameters for the standard case, and we choose 
the physical properties as stated by Li et al. [29] (shale 
gas), and by de Souza [12], which simulated the flow 
of natural gas, using cylindrical coordinates, to assess 
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the pressure of a producer well (vertical). As a sim¬ 
plification, we took C " c +1 as constant and equal to C sc 
(we intend to modify this in future work). 

Table 2: Parameters for the standard case. 


Parameter 

Value 

Unit 

C(f) 

1.0 io- 6 

psr 1 

Csc 

0.7 

scf/psi 

k r and k z 

4.0 10~ 6 

Darcy 

L-p 

1,250 

ft 

L z 

40.0 

ft 

L w f 

40.0 

ft 

n r 

40 

- 

n z 

3 

- 

p ini and p° 

4,500 

psi 

Psc 

14.65 

psi 

Qsc 

-1.0 10 4 

scf/day 

R 

10.73 ft 3 psi/R Ibm-mol 

t"max 

375 

day 

tol 

1.0 10" 6 

psi 

T 

609.67 

R 

T 

-*- sc 

519.67 

R 

7 

0.6 

- 

<5 A t 

1.1 

- 


0.0001 

day 

^tmax 

10.0 

day 

T 

1.41 

- 

</> and 

0.12 

- 


4.1. Numerical verification 

We employed different meshes in the study of nu¬ 
merical convergence: Meshes 1, 2, 3, and 4 with n r = 
10, 20, 40, and 80 cells, respectively. We kept, in all 
simulations, n z constant and equal to 3. 

Figures 4 and 5 bring the results for the four 
meshes mentioned in the last paragraph. The results 
show the pressure in the well as a function of time, 
obtained considering a production time equal to 375 
days. For this purpose, we use specialized and di¬ 
agnostic graphics, respectively. By the way, in this 
work, on the diagnostic plots, continuous lines repre¬ 
sent pressure drop, and dashed lines represent Bour- 
det derivative [8]. 

From the figures, we can see that we achieve 
numerical convergence as the number of cells n r 
increases, with the consequent overlap of pressure 
curves. Then, as a consequence of the results ob¬ 


tained, the 40 cells mesh was adopted as the stan¬ 
dard (lowest computational cost). From the pressure 
curves in Fig. 4, we realize that we capture two typi¬ 
cal regimes: a nearly horizontal line that corresponds 
to the storage in the well and an inclined straight line 
related to the transient flow regime. However, the ef¬ 
fects of external borders are still absent, reflecting the 
fact that the pressure at the reservoir frontier (r = r e ) 
remains equal to the pressure p ini . From the Bour- 
det derivative (Fig. 5), we can more clearly distinguish 
the two distinct flow regimes. First, we observe the 
inclined line related to the wellbore storage and af¬ 
ter a straight line corresponding to the transient flow 
regime, without boundary effects. 



Fig. 4: Numerical convergence under grid 
refinement, specialized plot. 



Fig. 5: Numerical convergence under grid 
refinement, diagnostic plot. 

We also carried out tests to observe the behav¬ 
ior of the solution concerning the variation of the time 
step A t n+1 = <5 A f Af ra . We consider two situations for 
prescribed A t ini and A t max \ numerical simulations 


www.ijaers.com 


Page | 231 





















International Journal of Advanced Engineering Research and Science (IJAERS) 
h ftps ://dx. doi. org/10.22161/ijaers. 76.29 


[Vo 1-7, Issue-6, Jun-2020] 
ISSN: 2349-6495(P) \ 2456-1908(0) 


with 6At = 1, Figs. 6 and 7, and pressure curves cal¬ 
culated with a variable S A t {8At 7 ^ 1). Figs. 8 and 9, 
respectively. 



Fig. 6: Wellbore pressure variation as a function of 
time step, specialized plot. 



Fig. 7: Wellbore pressure variation as a function of 
time step, diagnostic plot. 



Fig. 8: Wellbore pressure variation due to the growth 
of At, specialized plot. 



Fig. 9: Well pressure variation due to the growth 
of At, diagnostic plot. 


As we can see from these figures, the results did 
not show significant variations concerning changes in 
the well pressure and pressure derivative when we 
employ these A t ini and St values. Based on this veri¬ 
fication, we established the standard values for the ini¬ 
tial time step, the growth rate, and the maximum time 
step of Table 2. As already mentioned, the use of a 
small initial time step contributes to enabling the anal¬ 
ysis of the pressure variation in the initial moments of 
production (typical in Well Test Analysis using numer¬ 
ical simulation). 

We also studied the effect of the tolerance used 
to stipulate convergence in the internal (Gauss-Seidel 
method) and external (Picard method) iterations. Fig¬ 
ures 10 and 11 present the results, which corroborate 
the choice of tolerance for the standard case since we 
did not detect significant differences in the results. 



Fig. 10: Wellbore pressure variation as a function of 
tolerance, specialized plot. 
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Fig. 11: Wellbore pressure variation as a function of 
tolerance, diagnostic plot. 


As an additional check, Figs. 12 (specialized plot) 
and 13 (diagnostic plot) show our results obtained 
with higher values for porosity and permeability, 0.2 
and 1.0 10~ 3 Darcy respectively, and without taking 
into account the phenomena of slippage, formation 
damage, and wellbore storage. 



Fig. 12: Well pressure variation for higher 

permeability and porosity, and Darcy’s flow, 
specialized plot. 


This simulation is the same as those performed by 
de Souza [12] for the flow considering only the clas¬ 
sic Darcy law, whose results were verified through a 
direct confrontation with those obtained with the com¬ 
mercial simulator IMEX/CMG [11]. IMEX is a com¬ 
mercial simulator widely known and used in reservoir 
simulation. 



icr 4 nr 3 icr 2 nr 1 io° io 1 io 2 io 3 

t (days) 


Fig. 13: Wellbore pressure variation for higher 

permeability and porosity, and Darcy’s flow, 
diagnostic plot. 

The values here numerically calculated accurately 
reproduce those determined by de Souza [12], There¬ 
fore, indirectly, it was possible to validate our simulator 
in the case of flow governed by the classic Darcy law, 
without the effects of slippage, formation damage, and 
wellbore storage. 

4.2. Sensitivity analysis 

We begin by the sensitivity analysis of the effects 
on the wellbore pressure when we change the reser¬ 
voir’s permeability. In Figs. 14 and 15, we can observe 
the specialized and diagnostic plots for wellbore pres¬ 
sure variation. 



t (days) 

Fig. 14: Wellbore pressure variation as a function of 
absolute permeability, specialized plot. 

The higher the permeability, the lower the pres¬ 
sure drop, following the modified Darcy’s law. It is 
worth mentioning that permeability has a significant 
impact on flow in porous media, and is often respon¬ 
sible for defining the economic viability of the reser- 
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voir, including decision making on a hydraulic fractur¬ 
ing operation. In the case of reservoirs, where the 
slip flow regime occurs, the permeability value also 
influences Kn, and the Knudsen number increases if 
we decrease it, with a consequent increase in the ap¬ 
parent permeability value. We also noticed that the 
higher the permeability, the shorter the duration of the 
transition for the transient regime. 



t (days) 

Fig. 15: Wellbore pressure variation as a function of 
absolute permeability, diagnostic plot. 

In this test, and for all curves, we note that bound¬ 
ary effects are absent. We can verify this by inspect¬ 
ing the end of the transient flow regime curve, in the 
specialized plot, and the curve corresponding to the 
transient flow regime in the diagnostic plot (Bourdet 
derivative). Besides that, we can highlight that well¬ 
bore storage physically depends on permeability in a 
way that for lower permeability, wellbore storage effect 
is more prolonged. 

Just as we did with permeability, we used different 
porosity values to check its influence in the wellbore 
pressure variation. We show the results in Figs. 16 
and 17. 

However, for the porosity values that we retained, 
we only note small variations in the well pressure 
curves (p W f). From the figures, we observe that 
higher porosity values lead to lower pressure drop 
for maximum production time, and this is due to the 
higher volume of gas that we can produce. There¬ 
fore, in this situation, the well pressure drops less for 
a fixed production flow rate. Porosity also influences 
the value of the Knudsen number but acting in the 
opposite direction. Regarding the boundary effects, 
the smaller the porosity, the shorter the transient flow 
regime. Moreover, as expected, more significant well¬ 


bore storage effects occur for lower porosity values, 
although we only capture slight differences for the val¬ 
ues used in simulations. 



Fig. 16: Wellbore pressure variation as a function of 
porosity specialized plot. 



Fig. 17: Wellbore pressure variation as a function of 
porosity diagnostic plot. 


On the other hand, changing the tortuosity leads 
only to slight changes in the pressure values for the 
range that we considered (see Figs. 18 and 19). 

We note from the results in Fig. 18 a higher pres¬ 
sure drop as the value of tortuosity increases. This is 
consistent with the fact that the higher the tortuosity of 
the medium, the lower the apparent permeability must 
be and, therefore, the higher the flow resistance. 
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Fig. 18: Wellbore pressure variation as a function of 
tortuosity, specialized plot. 



Fig. 20: Wellbore pressure variation as a function of 
storage, specialized plot. 


From the results shown in these figures, we real¬ 
ize that we did not detect significant differences in the 
Bourdet derivative for this input simulation data. 



Fig. 19: Wellbore pressure variation as a function of 
tortuosity, diagnostic plot. 

4.3. Well Test Analysis 

We now focus on the specific cases of well test 
analysis. The initial transient response, during a well 
test, comes from the pressure difference applied to 
the volume of fluid initially present in the well. 

Figure 20 shows the storage effect on pressure 
behavior in the well. The higher the storage effect, 
the smaller the representative range of the transient 
flow regime. This means that a well test must be long 
enough for the storage effects to disappear so that 
the flow recorded in the well test identifies the tran¬ 
sient regime of reservoir production. In practice, it is 
from the well test that we obtain the data for the char¬ 
acterization of the well-reservoir system. 


The diagnostic plot (Fig. 21) shows that the well¬ 
bore storage strongly influences the results in the 
initial times, in both the pressure drop and pres¬ 
sure derivative. The higher the wellbore storage, the 
smaller the pressure drop and the later the transient 
flow regime will be. 



t (days) 

Fig. 21: Wellbore pressure variation as a function of 
storage, diagnostic plot. 

We also notice significant results when compar¬ 
ing the pressure variation for the flow without (k Q = 
k) and with the slippage effects (k Q = f(Kn)k). 
The same simulation conditions and parameters were 
maintained for both models, varying only the property 
of interest. 

Initially, only the responses resulting from the slip 
flow effect are studied, Fig. 22. Despite presenting 
a similar initial behavior (storage and transition to the 
transient regime), the model taking into account the 
phenomenon of slippage presents a lower pressure 
drop as production time increases. We can explain 
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this by the increase in apparent permeability due to 
the change in the Knudsen number. 



Fig. 22: Well pressure variation depending on gas 
slippage, specialized plot. 


Based on the results in Fig. 23, we note that the 
use of the classical Darcy’s law (lower apparent per¬ 
meability) leads to wellbore storage more significant 
in magnitude and duration. 



Fig. 23: Well pressure variation depending on gas 
slippage, diagnostic plot. 


We present now the results corresponding to three 
different flow rates in Fig. 24. We only change the flow 
rate values: -5,000, -10,000, and -20,000 scf/day. We 
call attention to the fact that we show the results in the 
form of the pressure variation (A p w f) divided by the 
production flow rate (Q sc ), a graph applied in the area 
of well test analysis. We perceive non-linear behav¬ 
ior through the variation of pressure curves as time 
progress. 


•icr 2 



Fig. 24: Variation of the ratio Ap w f/Q sc as a function 
of flow rate. 


Figures 25 and 26 display pressure variation 
curves in a region damaged by contact with fluids ex¬ 
ternal to the reservoir during drilling and completion, 
and, as a consequence, we have a decreased per¬ 
meability in this region (it was applied the Hawkins’ 
formation damage model [39]). We assumed that the 
damaged region measures 2.77 ft in all tests run. 

After the beginning of the transient regime, we can 
verify a significant variation in pressure drop as the 
damage increases (the relationship between the per¬ 
meability of the damaged region, k s , and the non- 
damaged region, k). Furthermore, we observe that 
the effect of storage has a longer duration and magni¬ 
tude for higher formation damage values (as a result 
of reduced permeability around the wellbore). 



Fig. 25: Well pressure variation due to damage to 
formation, specialized plot. 
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t (days) 


Fig. 26: Wellbore pressure variation due to damage 
to formation, diagnostic plot. 

Finally, to show the appearance of the border ef¬ 
fects, even for a shale gas reservoir, we performed 
a numerical simulation using the data presented 
in Table 2, except for the following three values: 
tmax= 4,280 days, Q sc =-5,000 scf/day and L r =625 ft. 

Therefore, at the end of 12 years of production and 
after the transient regime, we realize that the pressure 
curve changes its inclination (downward curvature), 
Fig. 27, as a result of the emergence of border ef¬ 
fects. As we can see, the border effects can (in some 
cases) appear only after a long period of time. This 
fact shows the practical difficulty of carrying out well 
pressure tests, which can take from a few days, in 
general, to months when long-lasting. Finally, we can 
also detect the same effect in Fig. 28. It corresponds 
to the change in the curve slope in the final simulation 
times, after the horizontal line for the Bourdet deriva¬ 
tive. 



t (days) 


Fig. 27: Well pressure variation due to the 

appearance of border effects, specialized 
plot. 



t (days) 

Fig. 28: Wellbore pressure variation due to the 

appearance of border effects, diagnostic plot. 

In all tests performed, the Knudsen number kept 
its value at Kn < 10, outside the range of molecular 
flow. 

V. CONCLUSION 

We carried out a study to better understand single¬ 
phase flow in petroleum reservoirs, aiming to maxi¬ 
mize hydrocarbon recovery. The numerical simulation 
allows testing different production scenarios in less 
time so that we can choose an optimized production 
plan that leads to economically viable exploration. 

The results of this work showed the importance of 
applying a complete model (including the effects of 
slippage, formation damage and wellbore storage) to 
study single-phase flow in shale gas reservoirs. The 
use of the classic Darcy’s law can lead to well pres¬ 
sure results that do not correspond to reality due to 
the non-consideration of these effects. This is essen¬ 
tial information that must be taken into account if we 
are to have reservoirs producing and generating prof¬ 
its. 

From pressure tests in vertical wells, we could 
obtain the reservoir properties by solving an inverse 
problem using the physical model of this work. There¬ 
fore, depending on the results, the use of horizontal 
wells or hydraulic fracturing could make production vi¬ 
able for some reservoirs. 

As expected, the nonlinear behavior of the results 
was adequately detected, especially in cases where 
we varied the permeability values, which appear ex¬ 
plicitly in the governing partial differential equation. 
Also, we captured the wellbore storage and formation 
damage effects, expanding the scope of this study, 
and we verified the influence of these phenomena on 
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the appearance of the transient flow regime and bor¬ 
der effects. 
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